Conceptualization of seeded region growing by 
pixels aggregation. Part 4: Simple, generic and 
robust extraction of grains in granular materials 
obtained by X-ray tomography 
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Abstract — This paper proposes a simple, generic and robust 
method to extract the grains from experimental tridimensionnal 
images of granular materials obtained by X-ray tomography. 
This extraction has two steps: segmentation and splitting. For 
the segmentation step, if there is a sufficient contrast between 
the different components, a classical threshold procedure followed 
by a succession of morphological filters can be applied. If not, 
and if the boundary needs to be localized precisely, a watershed 
transformation controlled by labels is applied. The basement 
of this transformation is to localize a label included in the 
component and another label in the component complementary. 
A "soft" threshold following by an opening is applied on the initial 
image to localize a label in a component. For any segmentation 
procedure, the visualisation shows a problem: some groups of two 
grains, close one to each other, become connected. So if a classical 
cluster procedure is applied on the segmented binary image, these 
numerical connected grains are considered as a single grain. To 
overcome this problem, we applied a procedure introduced by 
L. Vincent in 1993. 

This grains extraction is tested for various complexes porous 
media and granular material, to predict various properties 
(diffusion, electrical conductivity, deformation field) in a good 
agreement with experiment data. 

Index Terms — dynamic, mathematical morphology, segmenta- 
tion, watershed. 



I. Introduction 

The 3D reconstruction of microstructures by microtomogra- 
phy is based on X-ray radiography. From several radiographs, 
obtained after rotations of a specimen, and with an appropriate 
algorithm, 3D images of the sample are obtained at a micro 
scale. This article opens the way to predict the physical 
properties using 3D images where each voxel belongs to a 
grain of the granular phase or to the matrix phasqj. For 
example, 3D images of the microstructure can be introduced 
in a numerical homogenization process by computer, in order 
to predict their macroscopic physical properties from the 
microscopic properties. The first step to achieve this goal is 
to perform a correct grains extraction. Due to the acquisition 
of radiographs and to the numerical reconstruction, noise is 
apparent in the images, as well as some linear grey-level 
artefacts resulting in impressive images with generally a too 
weak quality for a quantitative and automatic use. 
This paper proposes a simple, generic and robust method 

! This matrix phase can be void (granular A and D) or solid (granular B 
and C). 



to achieve this goal. Simple means that this method can be 
used by anybody who is not a specialist of image processing. 
Generic means that this method can be applied in a wide 
range of materials. This method has been applied for granular 
materials (see figure [B but its extension to other materials 
is straightforward. Robust means that the extraction is few 
sensitive with a "little" variation of the parameters. 
This method has two steps: 

1) segmentation: to transform the grey-level image into an 
image with different components, 

2) splitting: to transform a granular component into a set 
of grains. 

1) Depending on the histogram shape, the segmentation is 
applied on the initial image to use the tint information or on 
the gradient image to use the boundary information as sharp 
variation of grey-level. In these both approaches, a succession 
of morphological filters is used: 

• after the threshold, to remove the small islands and to fill 
the holes of islands for the tint information, 

• before the watershed transformation, to localize two la- 
bels: one included in the component, the other included 
in the component complementary for the boundary infor- 
mation. 

2) Whatever the approach, the final result shows one difficulty: 
some groups of two grains close one to each other become 
connected after the segmentation. So if a classical cluster 
procedure is applied on the segmented binary image, these 
numerical connected grains are considered as a single grain. 
To overcome this problem, a splitting procedure is applied to 
separate the connected grains. 

The outline of the rest of the paper is as follows: in Sec. II, 
the images characteristics obtained by X-ray tomography and 
the different materials are explained. In Sec. Ill, the classical 
threshold procedure is presented. In Sec. IV, the watershed 
transformation controlled by labels and the localization labels 
are explained. In Sec. V, a splitting procedure is applied to 
individualize each grain in the granular component. In Sec. VI, 
concluding remarks are made. 



II. Materials and methods 

A. Microtomography 

Microtomography is a non-destructive 3D-characterisation 
technic providing a three-dimensional imaged. Each voxel of 
the image is associated to a cube included in the material, 
under investigation 1]. In first order, its grey-value is the 
space average of linear X-ray absorption coefficient of the 
different solids and fluids contained into it. But since more 
often the tomographic reconstruction amplifies the noise of 
the projections, and generates artefacts, there is extra- term 
given impressive images with generally a too weak quality 
for a quantitative and automatic use. Also, the materials are 
different in the chemical composition and in the geometrical 
organisation (see figure [T]). Due to the materials variety and 
the images defects, a generic, simple and robust segmentation 
procedure has been developed. 

B. Materials and applications 

For the granular A, the segmented data come from a 
mechanical triaxial test on a sand specimen realised under 
a synchrotron microtomograph (ESRF, ID 15 A) to follow the 
structural evolution of the granular media. Digital Image 
Correlation is used to observe and detect the strain localisation 
mechanisms at the grain scale[8], [12]. This work is funded 
by the French project ANR-05-BLAN-0192 (J. Desrues project 
coordinator). 

For the granular B and C, the segmented data come from py- 
rotechnical specimen realised by Laboratoire 3S-R, Grenoble 
under a laboratory microtomograph designed by Sky scan®. 
The finite element code is used to compute stresses and strains 
as well as other fields like the thermal flux and the temperature 
distribution in the material at a micro scale in using [11]. This 
work is funded by les Mines de Paris and le Centre d' Etudes 
de Gramat (DGA). 

For the granular D, the segmented data come from estail- 
lade limestone realised under a laboratory microtomograph 
designed by Skyscan® to understand the effects of the porous 
structure in resistivity index curves [9]. 

C. The image size 

The segmentation procedure must be efficient because the 
size of the images obtained by X-ray tomography is huge. 
The previous articles [20], [18], [19] enable to improve the 
efficiency of the algorithms used for this prupose. Typically, 
for the biggest image of this article, the granular D with a size 
equal to 700*700*700=0.35 Giga voxels, the segmentation 
requires less than 6 hours of executing time with an Intel(R) 
Xeon(R) CPU 3.00GH and the allocation of RAM is 16 GO. 
This reasonable time makes the practical use manageable on 
conventional computers. The algorithms are developed and 
implemented on the Open source software, called Population 
available at http://pmc.polytechnique.fr/~vta/Pop.tar 

2 Each image is a rectangular of the space E. E is a 3-dimensional discrete 
space Z n , consisting of lattice points which coordinates are all integers. The 
elements of a 3-dimensional image array are called voxels. 



D. Images characteristics 

For every image, the grey-level is coded on one byte (0- 
255) and a median filter has been applied to minimize the 
ring artefact and to smooth in keeping the sharpness of the 
boundary. The images have an uniform illumination and each 
component on the image has a specific brightnesfl For the 
visualization convenience, the results are sometimes presented 
in 2D but the method has been applied efficiently in 3D for 
all materials. 

III. Segmentation 

For each material, depending on the histogram shape, the 
classical threshold segmentation can be applied to extract a 
component, using tint information (1). If the contrast between 
the component and the background is low and if the boundary 
has to be well localized, the watershed transformation con- 
trolled by labels is applied using the boundary information 
(2). 

For the both approaches, a combination of morphological 
filters has to be applied in order to: 

1) match the visual segmentation for (1) (the combination 
is an opening followed by a closing), 

2) localize two labels for (2) (the combination is just an 
opening). 

This section is decomposed into two parts: threshold seg- 
mentation using tint information and watershed transformation 
using boundary information. 

A. Threshold segmentation using tint information 

1 ) Threshold operation: Given that each component has a 
specific brightness, the threshold operation uses this contrast 
information to extract the components. Threshold selection is 
usually based on the information contained in the grey level 
histogram of the image (see figure [2]). 

The correspondence between the component and a mode in 
the histogram is required for threshold segmentation. A certain 
threshold range is chosen. The label '1' is assigned to each 
image voxel, which gray-level belongs to this range, and '0' 
label is assigned to each image voxel which gray-level does 
not belong to this range. The threshold range is selected to best 
separate the mode of the histogram, see [16], [5], [14], [22]. 
Manually, the range is chosen on the valleys of the histogram 
(see figure [2]). For the granular A and B, the binary image 
seems like the real grains aside some holes and isolated voxels. 
The next paragraph explains how to correct these defects. For 
the granular C and D, the binary image is very different of 
the real grains. The reason is that one mode in the histogram 
is populated with voxels that form the two classes of grains 
because the contrast between these two classes is low. For 
these two granular materials, it is impossible to filter the binary 
image to get a good segmentation. The subsection IIII-BI gives 
a method to treat this task. 

3 This last assumption is not always verified. For example, the large grains 
with a medium average grey level in the granular B is divided into two 
components which chemical composition is different and which linear X-ray 
absorption coefficient is the same. Without more information, we consider 
these two components as one component. 



2 ) Morphological filtration: Briefly, the four basic operator 
of the mathematical morphology are presented[17]. 
Erosion of object A by the structural element B is defined 
b>0: 

AQ B = {z\(B) z c A} 

Dilatation of object A by the structural element B is defined 
by: 

A(BB = {z\(B) z HA^0} 

The Opening of A by B is obtained by the erosion of A by 
B, followed by dilatation of the resulting structure by B: 

AoB = (AQB)®B 

The Closing of A by B is obtained by the dilatation of A by 
B, followed by erosion of the resulting structure by B: 

AmB = (A®B)QB 

As the grains are isometric for all granular materials, the 
structural element is chosen isometric. The structural element 
is associated to the 26-connectivity in the cubic grid. 
To remove the isolated voxels for keeping the clusters size, the 
opening is applied. To fill the holes for keeping the clusters 
size, the closing is applied. Let A, the binary image after 
the threshold application. The filtration is only: (A o B) • 
B (see figure 0. We have a good agreement between the 
visual segmentation and the numerical segmentation but first 
the numerical boundary does not match closely the visual 
segmentation, second some groups of two grains close one 
to each other become connected after the segmentation. So, 
the following method has to be applied if 

1) for the purpose of the application, it is necessary to have 
a good match between the numerical boundary and the 
"real" boundary, 

2) the contrast between the components is low. 

B. Watershed transformation using boundary information 

1) Double labels watershed: An efficient segmentation 
procedure developed in mathematical morphology is the wa- 
tershed segmentation [3], usually implemented by a flooding 
process from labels. We recall these tools in this section, based 
on our application for the segmentation of the 3D images of 
materials. 

Minima watershed: Any greyscale image can be considered as 
a topographic surface and all boundaries as sharp variations of 
the grey level. When a gradient is applied to an image, bound- 
aries are enhanced. When the topographic surface obtained 
from the gradient is flooded from its minima, the waterfronts 
meet on watershed lines in 2D, and on watershed surfaces in 
3D. A partition of the investigated volume is obtained, where 
the catchments basins are separated by the watershed surfaces. 
However, in practice, this merging produces an important over- 
segmentation due to noise or local irregularities in the gradient 
image, generating a set of uncontrolled and unwanted markers. 
To avoid this problem coming from too many minima, the 
image, /, is usually filtered. It is a composition of vertical 

4 the inner boundary using for visualisation is denned as: dA = A\(AQB) 



(with respect to the grey level) and horizontal filters, in 
order to individualize each grain with a single marker. This 
individualisation step is complicated [20]. It is the reason why 
the labels controlled watershed is used[2]. 
Labels controlled watershed: The labels controlled watershed 
is similar to the minima watershed, beside a catchment basin 
is associated to a label (see figure [4]). 

Component by component : Since the materials are different in 
the chemical composition and in the geometrical organisation, 
it is difficult to extract each component in the same procedure. 
It is simpler to extract component by component. Starting 
from the simplest component to extract, we proceed to extract 
the next simplest component each time. Notice that the last 
component is easily extracted because it is the complementary 
of the other components addition. 
Except the last component, the extraction procedure is: 

1) to localize two labels: one included in the component 
and the other in the component complementary (the next 
paragraph is dedicated to this task), 

2) to apply the Deriche's operator [6] on the initial image 
to get the gradient image, 

3) to apply the watershed transformation controled by 
labels on this gradient image with these labels (see 
figure [5]) . 

The component is the catchment basin associated to the label 
included in it 

2) Localize a Label in a component: Two labels have 
to be localized: one included in the component and the 
other included in the component complementary. Since the 
component complementary can be composed by more than 
one component, this procedure has to be operated on each 
complementary component (see figure 0. The procedure uses 
the tint information like the threshold segmentation: a thresh- 
old followed by a succession of morphological filters. The 
major enhancement of this approach is: it is not necessary 
that the label localization match the component contrary 
to the threshold segmentation. It is just sufficient to have 
some labels in each connected component. The succession of 
morphological filter is just an opening. The opening operator 
has one parameter k such as 

A9 k B = (A@B k )GB k 
times k 

with B k = (B®B)...®B 

In this article, the selection of the threshold/opening param- 
eters and are done manually following these constraints (see 
table 0: 

1) the material specialist checks if the visual segmentation 
matches the numerical segmentation, 

2) if there is some experimental data about the volume 
fraction, we impose the correspondence between the 
experimental value and the numerical value obtained by 
segmentation. 

This manual limitation is attenuated by a good property: some 
small parameters modifications have no consequence on the 
final segmentation (see subsubsection IIII-B.31) . So it is easy 
to find the right parameters for a good segmentation because 
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the range of the right parameters is large. This simple method 
gives some good results for the four granular materials. The 
figure [6] shows the different steps for the extraction of one 
component for the granular A, B and C. The figure [7] shows 
the 3D visualization of the multi-component extraction. In 
the next sub sub section, a method to evaluate the robustness 
is presented, in more this method opens up the opportunity of 
an automatic evaluation of the parameters. 
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TABLE I 

B 1 AND B2 MEANS THE COMPONENT NUMBER 1 AND 2 OF THE 
COMPONENT COMPLEMENTARY. 



3) Robustness: A definition of the segmentation robustness 
can be the stability of some morphological descriptors depend- 
ing on the parameters fluctuation. This subsubsection has two 
paragraphs: morphological analysis and robustness evaluation. 

Morphological analysis: The /i-chord distribution 
function and the two-point correlation function provide a way 
to get a morphological analysis of the mesoscopic divided 
material and particularly for granular materials [4], [13], [15]. 
These functions can be computed either in 2D or 3D assuming 
some general conditions similar to the Hadwiger conditions. 
Their determination gives information about volume fraction, 
average granular size, mean curvature, granular shape, surface 
roughness and structural correlation. 

1) A chord is a segment belonging either to the grains or 
to the background and having its two extremities on the 
interface. As shown in Fig. [8l /i-chord are obtained by tracing 
random and homogeneous distributed straight lines through 
the microstructure. The chord-length distribution function 
gives the probability of getting a chord length between r 
et r + dr, belonging either to the grain, f 9 ^(r), or to the 
background f b ^(r). 

2) Let a yardstick drawn randomly in the material. The 
two-point correlation functions gives the probability of having 
both extremities of the yardstick belonging to the grains, 

4,2point( r ) or t0 the background f i 2 po int( r ) ( see 
figure O. 

Robustness evaluation: Let Sa(I) be the resulting binary 
image after the application of the threshold segmentation or 

5 For a biphasique material, the two point correlation functions of one phase 
is formulated in function of the two point correlation functions of the other 
phase. 



the double labels watershed on the initial image I with the 
parameters A. Let A5 the parameters that gives the "best" 
binary image. 

The method to check the stability of a segmentation method 
is: 

1) select a parameter, A, in the parameters, A, 

2) starting from well below of A5, increase step by step 
the chosen parameter, A^ = A;_i + e, until a value well 
above of A5, 

3) for each step, 

• apply the segmentation procedure to get the binary 
image Sa^I) associated to the parameter A^, 

• calculate the distance for the two point correlation 
function and for the chord-length distribution be- 
tween S Ai (J) and 5 Ai _ x (7)0. 

The figure [TOj) shows that the double labels watershed is 
more stable than the threshold segmentation. For the double 
labels watershed, plotting these curves will allow an automatic 
determination of the threshold value to localize the labels 
(see figure [TT]). Whatever the parameters, there is always a 
problem: some pairs of grains, close one to each other, become 
connected. A split procedure will be applied to separate the 
connected grains. 

IV. Split the grains 

For the both segmentation procedure, the final result shows 
one difficulty: some groups of two grains close one to each 
other become connected after the segmentation. So if a classi- 
cal cluster procedure [10] is applied on the segmented binary 
image, these numerical connected grains are considered as a 
single grain. To overcome this problem, we apply a procedure 
introduced by L. Vincent in 1993 [21]. 

The algorithm principle is to apply the watershed transforma- 
tion controlled by labels restricted by the binary image, 7, on 
the topographic surface that is the opposite of the distance 
function of /. In this approach, the crucial step is the labels 
localization. Each catchment basin is associated to a grain (see 
figure [T2l). 

In a first approach, the labels can be the regional minima of 
the opposite distance function of the binary image. However, 
in practice, this merging produces an important over-extraction 
of grains due to local irregularities of the binary image shape, 
generating a set of uncontrolled and unwanted labels. To 
avoid this problem coming from too many labels, the opposite 
distance function f of the binary image, 0(dist(I)), is filtered. 
A vertical filter is used in order to individualize each grain with 
a single label. The vertical filter uses the notion of dynamic 
with a parameter h. The operator is a geodesic erosion of the 
image 0(dist(I)) + h under the image 0(dist(I)) that fills 
the valleys with depth lower than h[7], [19] (see figure \\3\i. 
The splitting procedure has been successfully applied to the 
granular A-B-C0 (see figure [14]) because each grain of the 
granular phase has been correctly extracted. This method 

°The distance is: f J ■ . — f . . 2 
11 J g,2 point or n J g,2 point or /x" 

7 For the granular D, the physical property under investigation is the 

electrical resistivity and the segmented components image are sufficient to 

predict it. 
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has been applied for granular materials but its extension to 
other materials is straightforward. The figure Q3] shows the 
correspondance between the 3d grey-level image and the 
segmented image for the material presented in this paper and 
also for a cement paste. 

V. Conclusion 

X-ray tomography images involve specific problems of 
image segmentation, due to some artefacts generated during 
the reconstruction. In the case of granular materials, we have 
seen that many problems can appear during the segmentation. 
In the first part of this paper, we presented the classical 
threshold segmentation. To overcome the problem due to the 
low brightness between the components, we introduced a 
simple, generic, robust method of segmentation, called double 
labels watershed. This method is simple because the chain of 
operator is always the same, generic because this method has 
been applied in a wide range of material, robust because some 
"little" fluctuation of the parameters does not change the final 
segmentation. Then to split the connected grains generated by 
the segmentation in tridimensionnel image, we have applied a 
classical procedure introduced by L. Vincent. 
This simple, robust and generic segmentation can be applied 
efficiently in using the framework introduced in the previous 
articles. The initial goal, to give a method for any user not 
especially specialist in image processing, is reach. 
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Fig. 1. Figure 1: Granular Material A, size=450x450x200, resolution=14 
microns; figure 2: Granular Material B, size=500*500*100, resolution=3 
microns; figure 3: Granular Material C, size: 350x350x1 50, resolution=3 
microns, there are three components: the matrix in light and two grain classes: 
small dark grains and large grains with a medium average grey level ; figure 4: 
granular material D, size=700x700x700 resolution=3 microns, there are three 
components: the void in dark and two grain classes: one lighter and the other 
with a medium average grey level. 
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Fig. 2. Each serie is associated to a granular, the first image is the histogram, 
the second figure is the binary image after threshold, the third figure is 
the foreground of the binary image inner boundary on the initial image. 
Serie 1): granular A, the range of the thresold is 125:255, despite some holes 
in the numerical grains and some voxels outside the grains, the numerical 
segmentation is quite similar to the visual segmentation, serie 2): granular B, 
the range of the thresold is 0:125. Serie 3): granular C, the right mode is 
populated with voxels that form the two grains classes within the image. We 
have selected arbitrarily the range 120:145 to extract the darkest grains. The 
binary image shows some clusters of voxels corresponding to the middle grey 
grains, some lines corresponding to the boundary of the porosity (the darkest 
phase) due to a halo artefact and some isolated voxels corresponding to the 
light grains. Serie 4): granular D, the range is 80:144, idem that granular C 




Fig. 3. For both serie, the binary image is the application of threshold, the 
second image is the application of opening on the first image to remove the 
isolated voxels, the third image is the application of closing on the second 
image to fill the holes, last image: foreground of the inner boundary of the 
third image in red on the initial image. 
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Same color: mergin 
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Fig. 4. This transformation requires two images: the topographic surface 
(a grey-level image) and the label image. The process is: 1) association of 
each label to a hole 2) immersion 3) the water enters in the topographic by 
the holes and the catchment basins take the colour of the hole, 4) a part of 
the topographic is not merged although its level is under the level of the 
immersion, 5) fusion of two catchment basins of same colour, 6) creation of 
dam when two catchment basins have different colours. A video is available 
at http://pmc.polytechnique.fr/r^vta/water. mpeg 




Fig. 5. Starting from the initial image, two labels are localized: one included 
in the component and the other in the component complementary. To localize 
a label in a component, a threshold followed by opening is operated. Since 
the complementary is composed by n components (2 in this example), the 
localization is done for each component such as the final localization for 
the complementary is the addition of these specific localizations. When a 
watershed transformation is applied on the gradient with these two labels, 
the catchment basins are well localized on the grains. There is a good match 
between the visual segmentation and the numerical segmentation. 
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Fig. 8. A chord trough a granular material composed by two components 



Fig. 6. Visualization of each step for the granular material A, B and C: 
the first image is the initial image, the second is the label localization inside 
the component, the third image (and the fourth image for the granular C 
because the complementary is composed by two components) is the label 
localization inside each component of the complementary, the fourth image 
is the visualization of the double labels, the fifth image is the gradient image, 
the sixth image is the catchment basin associated to the label included in the 
component, the seventh image is the foreground of the inner boundary of the 
previous image on the initial image. 
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Fig. 7. The first image shows the white grains for the granular A. The 
second image shows the black grains for the granular B. The third image 
shows in yellow the small black grains and in blue the big grey grains for the 
granular C. The fourth image shows in blue the porosity (the black phase) 
and in yellow the black grains for the granular D. 




Fig. 9. A yardstick drawn randomly in a granular material composed by two 
components 
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stability of two point correlation function 
in function of threshold value 
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Fig. 10. Application for the granular material B: for the threshold seg- 
mentation, the threshold value is selected (the value 128 corresponds to the 
valley on the histogram) and for the double labels watershed, the threshold 
value to localize the label inside the grains is selected (the value 90 has been 
chosen manually to give a result matching the visual segmentation). For both 
distances, the double labels watershed is more stable of one decade than the 
threshold segmentation. 
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stability of chord lenght distribution 
in function of parameter fluctuation 
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mage after a dynamic filter 
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Fig. 13. The dynamic filter. Before the application of the dynamic filter, 
there are many regional minima (green bullets). After the application of the 
dynamic filter, there are only two regional minima. 
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Fig. 11. After a dilatation of parameter 5 on the graph, the minimum on 
both curves is 87. This automatic value is closed that the manual value: 90. 



I and three labels. 



D: the distance function of I 
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Fig. 14. Each colour represents a grain. For each serie, the first image is the 
visualization of the clusters on the segmented binary image and the second 
image is the catchment basins after the splitting procedure. The serie and the 
granular is 1-A, 2-B, 3-C. For each granular, the dynamic parameter is equal 
to 1. 



Fig. 12. Grains extraction: the watershed lines/surfaces are located on the 
narrow lines/surfaces because the narrow lines/surfaces are located on the 
crest lines of the opposite of the distance function of /. 
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Fig. 15. Corresdance between the 3d grey-level image and the segmented 
image. The last serie is a cement paste material composed by three phases: 
porosity (in black), hydrate (in middle grey) and anhydrous grains (in light). 
In the segmented image, the porosity is in blue and the the anhydrous grains 
is in yellow. 



